rm(list = ls())

library(tidyverse)
library(pbapply)

set.seed(02138)

## Get data

dat <- read_rds("data/Data_Clean.rds") %>% 
  mutate(year = as.numeric(year)) %>%
  filter(year %in% c('2017', '2021')) %>% 
  mutate(unit = coalesce(ags, vb_key))

## Get Treatment by AGS

treat <- dat %>%
  dplyr::select(one_of('unit', 'treat_categ')) %>%
  distinct(., .keep_all = T)

## Define outcomes 

outcomes <- c('cdu_csu', 'fdp', 'spd', 'afd', 'greens', 'left')

## Calculate first differences for all outcomes 

dat_fd <- dat %>%
  pivot_wider(id_cols = c('unit', 'land'),
              values_from = all_of(outcomes),
              names_from = 'year') %>%
  mutate(fd_cdu_csu = cdu_csu_2021 - cdu_csu_2017,
         fd_fdp = fdp_2021 - fdp_2017,
         fd_spd = spd_2021 - spd_2017,
         fd_afd = afd_2021 - afd_2017,
         fd_greens = greens_2021 - greens_2017,
         fd_left = left_2021 - left_2017) %>%
  dplyr::select(matches('fd_|unit')) %>%
  dplyr::select(-one_of('afd_2021', 'afd_2017')) %>%
  left_join(treat) %>%
  mutate(county_id = substr(unit, 1, 5))

## Combine schwer / sehr schwer treatment categories 

dat_fd <- dat_fd %>%
  mutate(treat_categ = recode(treat_categ,
                              'sehr schwer' = 'schwer')) %>%
  filter(!is.na(treat_categ))

## Land 

dat_fd <- dat_fd %>% 
  mutate(land = ifelse(substr(county_id, 1, 2) == '05', 'NRW',
                       'RP'))

## Function to get FD by sample

get_fd <- function(df) {
  
  ## Calculate grouped means 
  ## Cluster bootstrap for SEs 
  
  n_rep <- 1:2000 
  clusters <- unique(df$county_id)
  n_clusters <- length(clusters)
  
  simdat <- pblapply(n_rep, function(i){
    
    ## Sample Kreise with replacement 
    
    counties_sampled <- sample(clusters, size = n_clusters, replace = T)
    
    ## Create bootstraped sample 
    
    dat_bootstrap <- lapply(counties_sampled, function(cid){
      
      df %>% 
        filter(county_id == cid)
      
    }) %>%
      reduce(bind_rows)
    
    ## Calculate test statistic by group 
    
    test_stat <- dat_bootstrap %>%
      group_by(treat_categ) %>% 
      summarise_at(vars(contains('fd_')), mean, na.rm = T)  
    
    return(test_stat)
    
  }) %>%
    reduce(bind_rows)
  
  ## Calculate mean and CI
  
  boot_means <- simdat %>%
    group_by(treat_categ) %>%
    summarise_at(vars(contains('fd')), mean) %>%
    pivot_longer(cols = contains('fd'),
                 names_to = 'party',
                 values_to = 'change')
  
  boot_low <- simdat %>%
    group_by(treat_categ) %>%
    summarise_at(vars(contains('fd')), quantile, probs = .025) %>%
    pivot_longer(cols = contains('fd'),
                 names_to = 'party',
                 values_to = 'conf.low')
  
  
  boot_high <- simdat %>%
    group_by(treat_categ) %>%
    summarise_at(vars(contains('fd')), quantile, probs = .975) %>%
    pivot_longer(cols = contains('fd'),
                 names_to = 'party',
                 values_to = 'conf.high')
  
  ## Combine 
  
  plot_dat <- boot_means %>%
    left_join(boot_low) %>%
    left_join(boot_high)
  
  ## Clean up for plotting 
  
  plot_dat <- plot_dat %>%
    mutate(outcome_label = recode(party,
                                  'fd_cdu_csu' = 'CDU/CSU',
                                  'fd_fdp' = 'FDP',
                                  'fd_spd' = 'SPD',
                                  'fd_afd' = 'AfD',
                                  'fd_greens' = 'Greens',
                                  'fd_left' = 'Left party'),
           treat_label = recode(treat_categ,
                                'none' = 'Unaffected',
                                'schwach' = 'Weakly affected',
                                'schwer' = 'Highly affected')) %>%
    mutate(treat_label = fct_relevel(treat_label,
                                     'Highly affected', 'Weakly affected')) %>%
    mutate(outcome_label = fct_reorder(outcome_label, change)) %>%
    mutate_at(vars(one_of('conf.low', 'conf.high', 'change')), ~ .*100)
  
  plot_dat
  
}

## Do by land

res_all <- get_fd(dat_fd) %>% 
  mutate(ss = 'Full sample')
res_nrw <- get_fd(dat_fd %>% filter(land == 'NRW')) %>% 
  mutate(ss = 'NRW')
res_rp <- get_fd(dat_fd %>% filter(land == 'RP')) %>% 
  mutate(ss = 'RP')

## Combine

res_combine <- res_all %>% 
  bind_rows(res_nrw) %>% 
  bind_rows(res_rp)

## Spell out state names

res_combine <- res_combine %>%
  mutate(ss = recode(ss,
                     'NRW' = 'North Rhine-\nWestphalia',
                     'RP' = 'Rhineland-\nPalatinate')) 


## Plot This 

pd <- position_dodge(.4)

p1 <- ggplot(res_combine %>%
               filter(outcome_label == 'Greens'), 
             aes(x = ss, 
                 y = change, 
                 ymin = conf.low, 
                 ymax = conf.high,
                 col = treat_label,
                 shape = treat_label)) + 
  geom_point(position = pd) + 
  geom_errorbar(width = 0, 
                position = pd) + 
  theme_bw() + 
  theme(legend.position = 'bottom') + 
  labs(x = '',
       y = 'Average change in Green party\nvote share in p.p. (2017 - 2021)',
       shape = 'Flood intensity',
       col = 'Flood intensity') + 
  scale_color_grey()

p1


